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Nomenclature 


= principal moment of inertia of body about i coor- 
dinate 1^ + 2mr 2 , I 2 , ■ Ij + 2mr^ , respectively 

= principal moment of inertia of rigid core body about 
i body - fixed coordinate 

= stiffness coefficient characterizing nonrotating 
boom stiffness 

= tip mass of boom 

= applied torque about i coordinate 

= time 

= differentiation with respect to real time t 

= skew symmetric mode of elastic deformations 

= displacement of m tip mass from spinning steady 
state in i direction (m = I,Il) 

= perturbation [about spinning steady state) velocity 
about i coordinate 

= steady-state boom dimension in 2 -axis direction 
from center of mass to tip mass 

= the asymmetry in the setting of the booms 

= ratios of inertia and (I^-I^)/^ 

respectively 



= ratio of inertia 


o = 


1 + K, 


1 - K- 



Y = Zmr^/I-^ 

A. = d./mfi 
1 1 

Pi - u./2r 2 

o? = Kj/mn 2 


Ms'ssI&T fia’J v ) • ■' 

_X = fit 


fi 

w. 



( ) = d/dx 


2 

1 


= dimensionless inertia ratio 


= dimensionless damping ratio 

= general skew symmetric coordinate 

= dimensionless natural frequency coefficient of 
boom 

= dimensionless length ratio 

*= dimensionless time 
“ steady-state spin rate about 3 axis 
= dimensionless wobble ratio (i = 1, 2, 3) 

= differentiation with respect to x 


subscript i 


y 3 2mr 2 // ' I 3 



+i 


= index referring to three body- fixed coordinates 
(i = 1, 2, 3) 


= dimensionless inertia ratio 


= dimensionless inertia ratio 



= angular rotations about i-tK coordinate axis 


v 
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1. INTRODUCTION 


This report is a direct continuation of the decomposition- aggregation sta- 
bility analysis of the spinning Skylab presented in reference [1] . Therefore, 
for the mathematical basis, more detailed derivations of the state model, and 
computer implementations of the decomposition-aggregation analysis, this report 
will rely heavily on Jl] . 

In the preceeding report [1] , the linear equations of the Skylab motion, 
which include both the passive stabilization by extendable booms with .tip masses 
and the active stabilization by control torques about the body fixed axes, were 
decomposed into two sets of equations describing the wobble motion and the spin 
motion. Then, two sets of equations were treated as subsystems which were inter- 
connected by a coupling parameter representing the asymmetry in the booms - setting. 
Stability properties of each subsystem. were aggregated into a single quadratic 
Liapunov function. The vector Liapunov function was formed which had subsystem 

Liapunov functions as components. A linear second order comparison system was 

■> 

constructed in terms of the vector Liapunov function. Stability conditions of 
the aggregate comparison system provided estimates of the coupling parameter. 

Such an investigation is motivated by the fact that the mathematical model of 
the system is of high dimension and a straightforward analysis would become 
bogged down in the welter of detail requiring an excessive computer storage and 
time to complete the investigation. The multi-level decomposition-aggregation 
approach [2, 3] offers to solve the stability problems "piece-by-piece" and not 
only make more economical the computer use, but also reduce the liability of the 
errors in the analysis. Furthermore, by decomposing the system into parts that 
have important physical meaning, the decomposition-aggregation approach yields 
significant structural information about the behavior of the system, which is 
not generally available in a straightforward stability investigation. 



The decomposition- aggregation stability analysis is based upon the 
Liapunov stability theory, and it is inherently conservative. Therefore, 
the stability region of the structural parameter obtained in [1] is rela- 
tively small. 

The purpose of this work is to considerably improve the estimates of 
the stability regions by formulating and resolving a proper maximization 
problem [4, 5]. The solution of the problem provides the best estimate of 
the maximal value of the structural parameter and at the same time -yields 
the optimum comparison system, which can be used to determine the degree of 
stability of the Skylab. The analysis procedure is completely computerized 
resulting in a flexible and powerful tool for stability considerations of 
large-scale linear as well as nonlinear systems. 

The research reported herein was performed by S. M. Cuk under the super- 
vision of D. D. Siljak, and was used as a part of S. M. Cuk's M.S.E.E. Thesis, 
at the Electrical Engineering and Computer Science Department, University of 
Santa Clara. 
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2. PASSIVE STABILIZATION 


In this section, we will consider the passive stabilization of the spinning 
Skylab by extendable booms attached to the body of the vehicle. The linear vec- 
tor state equation describing the vehicle 

S: x'cO = PxCO (2.1) 

is obtained in fl] from the linearized equations of motion 


wobble J 
motion 1 


spin J 
motion i 


I 1 W 1 + (*3 “ l z ) ^2 + mr 2 ( u 3 + Q u 3^ 

-mr 3 (2^ + ii 2 - U u 2 ) = 0 

( J 1 " *3) ^1 + *2 ^2 + mr 3 (^1 * “ ' 2fiu 2 ) = 0 

i # 2 

2mr 2 (w 1 + nw 2 ) + m ii 3 + d^ u 3 + (k^ + mq ) u^ « 0 

i 

i 

i 

i ! 

2mr^ (ftw^ + Wg) + mii^ 

+ d-^ u 1 + u x - 2mqu 2 = 0 
2mr^ (-w^ + nw 2 ) + 2mftu^ + mii 2 

+ d 2 ^2 + (^2 ' m ^) u 2 = 0 ’ 


( 2 . 2 ) 


(2.3) 


derived in [4]. The symbols in equations (2.2-3) are introduced in the Nomen- 
clature. 

An important feature of equations (2.2-3) is that when = 0 , they be- 
come uncoupled into two sets of equations: the wobble motion (w^, w 2> u^) des- 
cribed by (2.2); and spin motion (w^, u-^, u 2 ) described by (2.3). The influence 
of the asymmetry in the arrangements of the booms (r_ j- 0) will be treated as 
the coupling parameter between the two motions. In the decomposition-aggregation 
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analysis, each motion represents a subsystem, and the coupling parameter can be 
made to appear explicitly in the interconnections among the two subsystems. 
Passive control equations (2.2-3) can be rewritten as follows: 

f If f If 

V 1 ‘ K 1 v 2 + + P3) " + P 2 “ = 0 

t II I 

-K^av^^ + av 2 + Sy^ ~ v i ~ ^2) = 0 

V 1 + v 2 + P 3 + A 3 y 3 + ^ a 3 + 1) P3 = 0 ( 2 - 4 ) 

1 f t ! 2 1 

5(v, + v 2 ) + y, + + cr^y^ - 2y 2 = 0 

5('vj + v 2 ) + 2y^ + y 2 + & 2 y 2 + (0^ - 1) y 2 = 0 , (2.5) 

where the notation is again as in the Nomenclature. The dimensionless parameter 
5 = the coupling parameter between the two sets of equations (2.4) 

and (2.5). 

The state space representation (2.1) of the over-all system S described 
by (2.4-5), is obtained by choosing the state 8-vector x(x) as 

? f f T 

X ( T ) = (v^ \>2 y-^ y i v 2 v 2) ’ (2.6) 

The system S of equation (2.1) can be decomposed into two interconnected 
subsystems described by 

^1’ x i^ T ) 

~*2 ’ X 2 


wobble 

motion 


spin 

motion 


= P^Ct) + £ 2 Q n (S) x x (t) + ?Q 12 (0 x 2 (t) (2.7) 

= P 2 x 2 (t) + CQ 21 (€) x x (t) + 5 2 Q 22 (?) x 2 (t) (2.8) 
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where the state vectors x(x) , x^(x) , x^Cx) of the system S and the two sub- 
systems and S 2 are 


x(x) = 


x x (x) 

x 2 M 


X 1 W " Cv l v 2 

I 

x 2 (t) = (y.^ y-^ 


\T 

>*3 “3 1 


\T 

y 2 v 2 ) » 


(2.9) 


In (2.7-8), the 4x4 matrices and P 2 correspond to the "free" 

subsystems and S 2 , and the 4 x 4 matrices Q^(£) , Q 12 (£), Q 2 ^(C)> 
and Q 22 (C) represent the interconnections between the two subsystems. 

In order to extract the subsystem matrices P^ and P 2 independent of 
the coupling parameter £ and obtain the decomposition of (2.7-8), it was 
necessary to use the following identities : 


i___ = I + 

2 a 


Z 2 y 


a - £ y 
1 


a(a-£ y) 
1 


Y - 1 + 5 


£ 2 y 


Y - 1 


(1-Y) (y-Hy) 


( 2 . 10 ) 


The matrices of (2.7-8) are: 

Pl2 " - P42 = (K 1 + M 1 ‘ 

P 13 “ Y °3' ( 1 ' ' 77 

P 14 = A 3 y /( 1 " Y ^ 

1 _ „ ■ 

P 21 K 2 

p 43 = - (03 + 1 - y )/(1- y ) 

P 44 = - a 3 /(i - Y ) 


p l = 


0 

1 

P 21 

0 

0 


Az 

0 


0 

1 

p 42 


a3 

0 

0 

4 


■14 

0 

1 

& 


5 



P 2 


0 1 

2 2 
P 21 P 22 


0 

0 


0 

2 

P 42 


0 

0 

0 

P 43 


0 

2 

p 24 


44 


2 2 
P 21 " ' °1 

P 22 = - A 1 

P 24 “ 2 

P 42 = ‘ 2 

p 43 = ‘^2 ' ^ 

P 44 = ‘ A 2 


QnCO - 


n li li ill 

0 q 12 q 13 q 14 


ll 


q ii 0 0 0 


0 0 


Q 12 0O - 


n 

*12 


11 

q 13 


q 14 


„11 _ C2Y + K x - 1) 

q 4? ~ ? 

(1-Y)(1-Y-C Y) 


2 2 
Y a- 


= -q 


‘ 43 ci-y)(i-y-S 2 y) 


n 11 - 

' q 44 


y 2 a, 


(i-y>(i-y-5 y) 


0 

11 

q 42 

11 

q 43 

11 

q 44_ 

q 11 . 

(k 2 * 1) 


21 

a-? Z Y 

0 

12 

q 12 

_12 

q 13 

0 

/ 

12 - 
q 12 

-qj 2 = 2 y/(1-Y-? 2 y) 

12 

q 21 

12 

q 22 

0 

0 

12 

q 13 

" q 43 = -y/C 1_ Y'? 2 y) 

0 

0 

0 

0 

12 _ 
q 21 

yO 2 + D/Co - ? 2 y) 

0 

12 

q 12 

12 

q l3 

0 

12 - 
q 22 ~ 

YA-j/Co - c 2 y) 
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Q n ft) 



<*21 = _ C 1 + K 2 ) a /(“ " 

qj 2 = (2a + K x - 1)/(1 - y -5 2 y) 

21 2 /n 2 
q 43 Y< V 

q 44 = Y4 3 /(1 -y-5 2 y) 




22 , 2 A , w , _2 % 
q 21 = " Y ^ a l ’ 5 Y) 

q 21 = " YA /^ a " ^y) 

qj 2 = 2 Y /(1-Y-C 2 Y) 

-Y/d-Y-5 2 Y) 


( 2 . 11 ) 


The structural configuration of the system 5 as composed of the two 

subsystems 5^ and and the interconnections between them through the 

coupling parameter 5 can then be depicted as in Fig. 2.1a. 

It is obvious that the system of Fig. 2.1a becomes that of Fig. 2.1b 

when 5 = 0. When 5 -*■ °° , the system of Fig. 2.1a is again decoupled into 

the two subsystems shown in Fig. 2.1c, because the interconnection matrices 

2 2 

5Q 12 (5), 5Q 21 (5) and 5 Q-q( 5), 5 Q 22 (0 become zero and constant matrices, 
respectively. 
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The subsequent stability analysis shows that the free subsystems and 
S 2 (£ = 0) are stable. It is easy to check, however, that the decoupled sub- 
systems in Fig. 2.1c are unstable. Therefore, our main objective is to deter- 
mine the best estimate £^ of the maximum allowable value of £ which lies 
between the two extremes £ = 0 and £ = <*>, and for which the overall system 
of Fig. 2.1a is stable. 

On the basis of the Skylab physical characteristics the matrices Q^(£) 

(i, j =1, 2) , can be made independent of £ and denoted by . This is 

2 2 

accomplished by neglecting the term £ y = 0.197£ with respect to the terms 
1 - y = 0.803 and a = 5.52 since £ << 1 . 

After this simplification, the numbers £.. (i = 1, 2) of the norm of the 
coupling matrices . can be computed using 

- UCQy Qij j] 1/2 . i, j = 1, 2 (2.12) 

The subsystem Liapunov functions v^, are chosen as 

v i Cx) = (x^i) 172 , i = 1 > 2 (2-13) 


and the differential inequality 


v < Av 


(2.14) 


is formed following [1-3] , where v = (v^ v^) is the vector Liapunov function 
and the aggregate 2x2 matrix A 


A = 


x A(G 1 ) 


2 AOLjT 


+ £ 


2 A 0V 
HH^T 


(a^j) is defined as 

g 12 [g| aCH 2 ) 
^ 1/2 (H 1 ) A 1/2 (H 2 ) 


?2l ! ^ I A (^2^ 
7^ i 1/2 (H 2 ) 


A (GO 


2 a TH jr + 5 xiry 


2 A ^ H 2^ 


(2.15) 
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In (2.15), A and A denote minimum and maximum eigenvalues of the in- 
dicated matrices, respectively. The subsystem Liapunov functions in (2.13) 
are obtained by solving the Liapunov matrix equation 

P T H. + H-P. = -G. , i = 1, 2 (2.16) 

using the direct method of solution. From [2, 3], the conditions 

a>^ < 0 , det A > 0 , (2.17) 

are necessary and sufficient for stability of A in (2.15), and sufficient for 
stability of the overall system (2.1). For every choice of positive definite 
symmetric matrices G^, G^ in (2.16) inequalities (2.17) yield an estimate of 
the stability region of the coupling parameter £ . The special choice of G^ 

= G^ = I in (2.16) produces positive definite matrices H^, H 2 and establishes 
stability of the decoupled subsystems 5^, 5 2 when £ = 0 [1]. 

The range of the coupling parameter £ obtained in [1] is small due to the 
conservativeness of the stability procedure. However, the maximum estimate of 
£ could be considerably increased by a proper choice of the matrices G^, i = 

1, 2, in (2.16). A meaningful optimization problem can be formulated as the 
maximization of £ over all matrices G^, G 2 . In Section 4, the corresponding 
optimization problem is formulated and resolved. [4,5] 
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3. ACTIVE STABILIZATION 


In order to inertially fix the axis of the Skylab pointed at the sun in 
presence of disturbance torques, attitude control torques must be applied to 
the vehicle [4] . The control torques depend on error signals that are propor- 
tional to the angle between the inertially fixed axis and the solar vector. 

Sun sensors and rate gyros on the present Skylab can readily provide the sig- 
nals <f>, , <f> 2 > w^ and needed for control. 

Again, the linear vector state equation is considered 

S: x'(t) = Px(x) , (3.1) 

which is obtained in [1] from the linearized equations of motion 


r 


wobble < 
motion 


I^i + (I 3 - I 2 )Ow 2 + mr 2^3 + n u 3-^ 

-mr 3 (2fiu^ + ii 2 - n u^) = T^ 

(1^ - I 3 ) nw n + I 2 w 2 + mr 3 (ii.^ - n 2 u.^ - 2ftu 2 ) = T 2 

Zmr^ (w^ + 0w 2 ) + mii 3 + d 3 u 3 + (k 3 + mfl Z ) u 3 =0 (3.2) 


spin / 
motion 1 


I 3 W 3 ’ mr 2 *- u l " Zn ^2' ) T 3 
2mr 3 (Ow^ + w 2 ) - 2mr 2 w 3 + mii^ 

+ d^ Uj + k^ u^ - 2mftu 2 = 0 

2mr 3 (-^ + ttw 2 ) - 4mr 2 «w 3 + 2mfJu 1 
2 

+ d 2 u 2 + (k 2 - mh ) u 2 = 0 , 


+ m(i 2 


(3.3) 
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derived in [4] using the linear control law 


T = a<j) + 8 w . (3-4) 

T T 

In (3.4), T = [T x T 2 T 3 ] is the vector of control torques; <j> = 

T 

is the vector of angular rotations; to = [w^ w^ + fi] is the vector of ang- 
ular velocities; a, $ are 3><3 matrices 



“ll 

a 12 

0 


e H 

*12 

0 

a = 

a 21 

a 22 

0 

, B = 

*21 

*22 

0 


_ 0 

0 

Oj 


0 

0 

*33 


and kinematic relationships are 



~1 

0 

0 “ 


~0 

fi 

0 ~ 

w = 

0 

-1 

0 

♦ + 

-fi 

0 

0 


_0 

0 

- 1 _ 


_0 

0 

0 _ 


The control law in this study is chosen. as 


(3.5) 


(3.6) 


a 


12 



all other a. . = 0 
il 


3 n = IjBfi 


B 33 = - I^fip, all other 



(3.7) 


so that the normalized control torques v = [v.. v ? v_] T = [T,/I,fi 2 T /I fi 2 

Jl Lt o i. Jl Z 1 

are 


v x = (e + 6 )<j> 2 - 6 ^ 



v 3 = p <t> 3 (3.8) 
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Referring to equations (3.2) and (3.3), the control torque T^ is used to sta- 
bilize the subsystem (wobble motion) , and the torque is used to sta- 
bilize the subsystem (spin motion) . 

In (3.8), e, 6, p are control parameters to be selected in the stabili- 
zation process. 

Upon introducing these transformations the linearized equations of motion 
become : 

-<t' 1 - ( 1 +K 1 )<f '2 * " yG-I-j+Vj) + ?y(2u 1 +y 2 -vi 2 ) + (e+6)<j> 2 - <54^ = 0 

(l + K-^)4>-^ + ^2 a< t>2 + 0l( i > 2 ' 5v(u 1 -y 1 -2y 2 ) = 0 

! 

‘ + y 3 + ^3^3 + (R3 + 1 )Uj = 0 ( 3 - 9 ) 

I 

It If I I 

B4> 3 + yG^ - 2p 2 ) + p4> 3 = 0 

i 1 

I i 

i it n it '2 1 

-2c<t> 1 - W 2 + ^3 + y i + Vl + °i y i ' 2y 2 + 5<f> 2 = 0 

5 ^ - 2 £< J >2 - 5* 1 + 2<> 3 + 2 y ^ + P 2 + A 2 y 2 + (o 2 " 1 ' )v 2 = 0 * ^ 3 ‘ 10 ^ 

The state space representation of the overall system S described by 
(3.9-10), is obtained by choosing the state 11-vector x(t) as 

• i i i i t t 

x(t) = (4^, <i>2» u 3 > <l> 1 , 4> 2 > v 3 , <l ) 3 > v 1 » Up y 2 » y 2^ ‘ C 3 - 11 ) 

The system S of equation (3.1) can be decomposed into two interconnected 

subsystems described by: 


wobble 

motion 


spin 

motion 
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wobb 1 p * 9 

motion V *lCt) = p x ^(x) + ♦ 5Q 12 CS)x 2 (x3 

spin S 2 : x 2 (t) = P, x,(t) + 5Q 71 C5)x-, (t) + C 2 Q„(5)x 7 ( j) 
motion ^ L 1 Ll L 


using the same procedure outlined in the previous section. 

The state vectors x(t) , (t) , x 2 (t) of the system S and two subsys- 

tems 5^ and are 


x(t) = 


—i ' ' ' x 

x i(t) x 1 (t) = <j> 2 y 3 ^ *2 y 3^ 


x 2 (t) x 2 (t) = 0 3 n x a 2 v 2 y 


In (3.12-13) the 6x 6 and 5x5 matrices P-^, P 2 correspond to the sub- 
systems 5^ and S 2 and 6x6, 6x5, 5x6 and 5x5 matrices (^(5), 
Q 12 (£), Q 2 2 (O > Q 2 2 (?) represent the interconnections between the two subsys- 
tems: 


0 

1 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

P 41 

1 

P 42 

1 

P 43 

1 

P 44 

1 

P 45 

0 

1 

P 52 

0 

P 54 

0 

1 

p 61 

1 

P 62 

1 

P 63 

P 64 

1 

P 65 
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P 41 = t K l + Y)/( 1- Y) 

P 42 = P 62 = “ (e + 6)/ (1 ~y) 
p 43 " " Yo 3 /(l-YJ 
P 44 = Pg 4 = «/(i-y) 

P« - P65 - ( M l)/d-v) 

P 46 * -v^/a-r) 


P 52 " K 2 

P 54 ' V 1 

Pgi = a+Kp/a-Y) 

P« + !-v)/a-T) 

P66 = -V»-^ 



2 

P 11 

2 

P 12 

4 

0 

0 



0 

0 

1 

0 

0 


2 

2 

P 31 

4 

2 

P33 

0 

2 



0 

0 

0 

0 

1 



-2 

0 

-2 

2 

P 54 

2 

P 55 


2 

P 11 

= - p/(0-Y) 


2 

P 32 

- 0oJ/(3-y) 

2 

p 12 

= ya^/Xe-y) 


P33 = 

- 3A 1 /(e-y) 

2 

P 13 

= yA x / (3-y) 


P54 = 

f 2 
- (a 2 

i- 1 ) 

2 

P 31 

= p/ (3-y) 


4- 

- A 2 
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12 12 
<41 = q 61 


12 _ 12 
q 44 q 64 


.. 2y 


lTiT 


2 

ya 2 


l-Y-5 Y 


12 

q 51 


y 3 p 


«(l"Y 3 )-€ 'y 


12 

q 52 " ~~z ~jr 

“Ci-YjJ-S Y 


12 12 YA 2 
q 45 = q 6S = 


l-Y-5 Y 


12 

q 53 


ya 


1 . 


°(!-Y 3 )-C Y 


Q 2 2^ " 


22 

l ll 


22 

q 31 


0 


g 22 

q 51 


22 

*12 


l 22 

l 32 

0 

0 


22 

*13 


22 

q 33 


0 

0 


0 


22 

q 54 


22 

q 55 


n 22 

^11 


PY: 


(1-Y 3 ) [a(l-Y 3 )-5 y] 


22 

q 33 


YAi 


(1-Y 3 ) [a(l-Y 3 )-5 Y] 


22 

q 12 


22 

q 13 


22 

q 31 


yy 5 ( ct i +1 -y 5 ) 

(1-Y 3 ) [«(1-Y 3 )-C 2 y] 
yy 3 A 1 ‘ 

(1-Y 3 ) [«(1-Y 3 )-5 2 y] 


PY 

3(i-y 3 ) [«(i-y 3 )-5 2 y] 


22 

l 51 


22 

q 54 


22 

q 55 = 


2y 


i-y-cy 


2 

Y°2 


(1-y-ey) 

r 

yA 7 


22 

q 32 


Y(a 1 +1-Y 3 ) 


(1-Y 3 )[o(1-y 3 )-5S] 
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n ,21 
0 q 12 

0 

21 

q 14 

0 

0 



0 0 

0 

0 

0 

0 


= 

n 21 

° q 32 

0 

21 

q 34 

0 

0 



0 0 

0 

0 

0 

0 



21 21 
q 51 q 52 

21 

q 53 

,21 

q 54 

21 

q S5 

21 

q 56 

21 

q 12 


r 3 a(K 2 + 1} 


a 21 = 

e+6 


(1-y 3 )-5 2 y 


q 52 

i-yVy 


21 

q 14 = 

y 2 

,(2«-l+K I ) 


a 21 = 

Y a 2 


a(l-Y 3 )-C Y 


q 53 

, f 2 

l-Y-5 Y 


21 

q 32 

«(y 1) 


n 21 - 

q 54 " 

5 


a(l-Y 3 )~5 2 Y 


1-Y-C" 

r~ 

y 

21 . 
q 34 

2a-K 1 -l 


a 21 = 

i-k^y 


«(1-y 3 )-5 2 y 


q 55 

i-yVy 


21 

q 51 ' 

1- 

\-2y 


a 21 = 

YA3 


1- 

J 
Y-5 y 


q 56 

1 -yVy 



The following identity relationships were used in order to get the sub- 
system matrices and independent of the coupling parameter 5 : 

1 _ 1 


cT 


ad-Yj 

1 1 


)-E, 2 C ^ 1 " T 3^ oc(1-y 3 ) [a(l-Y 3 )-5 2 y] 


6 2 y 


I-y-5 2 ’ 1 " y C1-y)(1-y-5 2 y) 


( 3 . 16 ) 
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The graphical interpretation of the interconnected subsystem and S 2 

is the same as in Fig. 2.1 of the previous section 

Again, on the basis of the physical characteristic of the Skylab given at 

the beginning of the report, the matrices Q—COfi, j = 1, 2) of (3.15) can 

be made independent of ? and denoted by Q— . This is obtained by neglecting 
2 -4 

the term £ y = 8.5 * 10 with respect to the terms 

1-Y = 0.803 and a(l- Y ^) = 5.33 . (3.17) 

Using the following specific values of the control parameters e, 6, p: 

e = 2.0 6 = -1.0 p = 1.0 (3.18) 

the same computational algorithm as in the Passive Control case can be applied 
to the stability analysis of the Active Control as shown in [1] . 

The choice of the 6x6 and 5x5 identity matrices for the G matri- 
ces of the first and second subsystem results in 6x6 and 5x5 positive 

i 

definite matrices H^, H 2 and establishes the global asymptotic stability of 
the subsystems. 

The internal of the coupling parameter £ obtained in [1] , for which the 

i i ' 

overall system is globally exponentially stable, is relatively small due to the 
following reasons : 

1. The inherent conservativeness of the stability analysis; 

2. The choice of the matrices G^, G^ is not the "best" 
regarding the maximum value of E, . 

As in the Passive Control case, the solution to the second problem (and, 
therefore reducing conservativeness of the analysis) lies in the solution to 
the corresponding maximization problem formulated and solved in the following 
Section 4. 
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4. MAXIMIZATION OF THE STABILITY REGION 


Now, several comments mentioned in previous sections about the conserva- 
tiveness of the obtained result can be clarified. First, up to this point, 
there has been no effort made to obtain the maximum possible value £ m of 5 . 
Secondly, there is a considerable flexibility in the choice of the matrices 
G^(i = 1, 2) which are only constrained to be symmetric and positive definite. 
This freedom in choice, however, if not used in an optimal way, might show as 
the big disadvantage. Furthermore, the following analysis will reveal the re- 
markable sensitivity of the size of stability region (that is £ ) on the choice 
of G^, G^ matrices. 

Therefore, the following question naturally arises: How to choose subsys- 
tem Liapunov functions such that the corresponding comparison stability in- 
equalities will not be overly sufficient. In other words, instead of having 
the comparison system as constructed in [1] , that is 

v x < (-0.96 x 10' 4 + 186.755 2 )v 1 + 1954.26 | C |v 2 
v 2 < 392.98|5|v 1 + (-11.54 x 10" 4 + 573.86c 2 )v 2 


-0.96 x 10' 4 + 186.75c 2 


1954. 26 | C 


A = 


(4.1) 


392.98 | ? | 


-11.54 x 10 4 + 573. 86C 2 


we would prefer the following comparison system 


20 



Vi 1 C-0J + e 11 5 ) V 1 + 0 12 I ^ l v 2 


v 2 - 0 21 


klv, + (~8' 


22 


5 2 )v 


2 


-«1 * e xV 2 


A 


o 





o 

6 22 C 



(4,2) 


in which the nonnegative numbers 0°, 0°^, 6^ > Q^i > 
, are chosen in such a way that 




dependent on , 


A° - A < 0 , VA G H 


where H is the class of matrices which satisfy conditions (2.17). 

Because v^, are Liapunov functions for the subsystems, v-^ >_ 0 , 
v 2 0 and the system (4.2) is then better than (4.1) as the comparison sys- 
tem since (4.2) implies (4.1) but the opposite is not true. Furthermore (4.2) 
gives better exponential estimate than (4.1). Obviously the system (4.2) will 
allow larger variation of £ for which the comparison system will be still 
stable. 

By introducing the norms of the coupling matrices and the absolute values 
of the coupling parameter £ we already "washed out" the fine structure of the 
interconnections, but with the bad choice of the subsystem Liapunov functions 
and therefore with the corresponding estimates n^, n^ 2 > \ 3 > (i = 1, 2) 
we could even worsen the situation. 

In order to enlarge the stability region with respect to the coupling 
parameter £ we can formulate the following optimization problem, that is 
- the problem of nonlinear mathematical programming: 
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PAobZm A. 


max |c | 

{G 1» G 2* 

T 

subject to the constraints: + = ~ G i ’ 1 = 1> 2 


VA G H (4.3) 


Here (Gp G 7 } is the set of all symmetric positive definite matrices 
for the two subsystems. Then for every choice of Gp Gp the corresponding 
Liapunov matrix equation should be solved for Hp H 7 and then the estimates 

n il’ n i2’ n i3’ n i4 can be computed from A((L), A0T), X(G i ), (i = 1, 2) [1] . 
Together with computed norms of the coupling matrices, these estimates would 
constitute the comparison matrix A. Checking the inequalities (2. 17) for dif- 
ferent values of £ would lead to a maximum E, for that particular choice 
of- Gp G^ . Since the optimization is carried out over {Gp G^} this prob- 
lem would be even more difficult than that of attacking the overall system 
directly by the same nonlinear programming approach. 

But following the same stream of ideas, we can reformulate the problem 

in the following way: If we neglect the higher order terms since the physical 

_2 

coupling parameter E, is assumed to be small of the order 10 we obtain the 
following comparison matrix A 


A = 


1 A(G X ) 

2 1(^1 


’ 21 ' 


A(f y. 

/xTHJT x(h 2 ) 



ACHjj ' 
'a(h J x(h 2 ) 

1 *( G 2 > 

2 ACHp" 


C4.4) 


This actually corresponds to retaining only cross -coupling terms in (2.7) and 
(2.8) and neglecting the self coupling terms as the second order effect. 
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Now, the first stability inequality a^ <0 is automatically satisfied 
and the second one gives 


2 i . UGj) X(H x ) X(G 2 ) A(H 2 ) 
g < T t : 2 1 

C 12 ? 21 A ^ A CH 2 ) 


(4.5) 


Since the norms \ ^ , 5 2 i t ^ ie intercoupling matrices Q l2 and Q 2 ^ 
are constant independent of the choice of matrices , G 2 the problem is 
then reduced to: 


Pno blm 8: 


subject to the constraints: 

Since the choices of G^ 
problem can be further decomposed into two problems of the same structure, 
that is: 

Problem C : 

1 A(G.) A(H ) 

max y * 

G. L A (H.) 

subject to the constraint: P T H. + H.P. = -G. , i = 1, 2 (4.7) 

1 x 1 1 i ’ * v ' 

The numerical solution of this problem may be still of the same kind of 


! V^) XCHj) a(G 2 ) a(h 2 ) 

max -r- -y~ ■ 

(G 1 , G 2 ) 4 A Z (H X ) A Z (H 2 ) 

P i H i + H i P i = -Gi , i = 1, 2 (4.6) 

and G 2 are independent of each other, the 
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difficulty as that of solving the problem directly without decomposition using 
mathematical programming techniques (for example, interior point-penalty func- 
tion approach). That is, to maximize £ subject to the constraint that the 
matrix P in equation (2.1) is stable. Therefore, nothing would be gained 
by decomposition-aggregation approach. The attempt to solve Problem C analy- 
tically involves, however, two rather serious difficulties: 

1. and IT are related through the Liapunov matrix equation 
(4.7) in a quite complicated manner. 

2. Dependence of the minimum and maximum eigenvalue of a symmetric 
positive definite matrix on its elements is numerical in nature. 

Despite these drawbacks, we will follow the second route because a better 
look at (4.7) in some special case of subsystems will suggest how to solve the 
problem in general. 

In order to illustrate how the conservativeness of the result can be 
strongly affected by the choice of the matrices G^(i = 1, 2) , the following 
simple example is analytically treated. 

Let us suppose that the subsystem matrices P^ and P^ to be given in- 
stead by (2.11) are 

P i = -I (i = 1, 2) (4.8) 

where I is, in this case, a 4x4 identity matrix; that would be the case 
when the subsystem states are completely decoupled. The Liapunov matrix equa- 
tion can then be directly solved in matrix form: 

-I T H. - H.I = -G. or 
ii i 

H i = I G i (i = l>2) (4.9) 
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From this equation the relationship between the eigenvalues of IT and 

G. is 
1 

X( H.) = i- X(G i ) 

A( H i )=i A (G.) (i = 1, 2) (4.10) 

Ptioblm C is then reduced to the following problem: 


max 


G. 

i 


A(G.) 

aTSTJ 


(i « i, 2) 


(4.11) 


This reveals how the final result very strongly depends on the ratio of 

minimum to maximum eigenvalue of the chosen matrix G^(i = 1, 2) . Actually 

having in mind original Pxoblm .8, dependence is to the fourth power of this 

ratio when G^ = G£ . For example, the choice of G^ = G 2 with X(G 1 )/A(G 1 ) 
-2 -8 

=10 would give 10 as the result. 

In this special case of subsystem matrices, expression (4.11) not only 
discovers the sensitivity of the problem, but also the explicit solution to 
the optimization problem. From (4.11) it is obvious that the maximum is ob- 
tained when X (G . ) /A (G . ) = 1 , (i = 1, 2) . This is the case if and only if 

g g 

G 1 = g l I ’ G 2 = g 2 l ’ *- g l’ g 2 > °- ) givin § H i = ^ 1 , h 2 = 1 and the 

maximum value 1 . Therefore, the identity matrices constitute the solution 
set to the PhobZem 8 in this special case of subsystems. 

But it would be misleading to conclude that identity matrix I as a 
choice for G^ (i = 1, 2) is the best in general as quite conservative 
results in reference [1] point out. Though this choice of matrix G^ is 
very common and the easiest one, it will be shown that it is far from an opti- 
mum one. 
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The preceding example, however, indicates the way to approach the solu- 
tion of quite formidable Pnoblrn C. Instead of looking at the problem in the 
original state space coordinates, we will first transform it to the domain of 
canonical coordinates, but in a special way. We will find the canonical trans- 
formations and for both subsystems and which will in case 

of complex eigenvalues + juf (i = 1, 2) of reduce it to the 

following canonical form: 



(4.12) 


or in the case of mixed real and complex eigenvalues of + j , 

°2 ’ °3 to t * ie canon i ca l form: 


0 0 


0 0 

0 0 a 2 0 


0 


0 0 



(4.13) 


The structure chosen for the canonical form induces a canonical basis of 
real vectors even when the eigenvalues of P^ and P^ are complex and there- 
fore T^ and T^ are always real. 

In (4.12) and (4.13) we have chosen 4x4 matrices only because of its 
relation to the Passive Control case and in order to simplify presentation. 
However, the next development and the solution of the optimization Problem C 
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will be carried out in general for an n x n matrix P-^ and q x q matrix 
P 2 . The only assumption which will be made is that either all the eigen- 
values of the subsystem matrices are of multiplicity one, or that the order 
of their multiplicity is the same as the number of linearly independent char- 
acteristic vectors corresponding to that multiple eigenvalue. This is, how- 
ever, quite a realistic assumption. 

Let us now choose the sign (tilda) to distinguish canonical basis, 
canonical forms, and transformed coupling matrices from the corresponding 
original one. 

By introducing the change of basis: 


or 


X 1 = T l x l “d x 2 = T 2*2 



1 

X 

1—* 

1 


1 

0 
rH 

H 

1 


1 

x < 

1— v 

1 


T x 0 

X = 


= 




= 



1 

X 

to 

1 


_° T 2 


1 

tvi 

'X 
1 


0 t 2 _ 


(4.14) 


(4.15) 


into the system description: 
wobble 


motion 

we obtain in canonical coordinates 


Sy x 1 P- L x 1 + E, Q 11 x 1 + ?Q 12 x 2 

(4.16) 

S 2 : X 2 P 2 X 2 + ^21 X 1 + 5 ^22 X 2 

(4.17) 


wobble . 

_ I 

motion 1 * 

X 1 

spin s . 

- { 
Y 

motion z' 

2 


r 2 P 2 T 2 X 2 + 5T 2 ^21 T 1 X 1 


’ T 2 ^22 T 2 X 2 


(4.19) 


or the following: 
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wobble 
motion ~1 


spin 
motion 


V 

*1 ^1 X 1 + ? ^11*1 + ^12 x 2 

(4.20) 

V 

X 2 = ^2*2 + ^21*1 + 5 ^22*2 

(4.21) 


Let us now show how the corresponding Liapunov matrix equations are 
changed under these transformations: 

Since 


P. = T^P.T. and pTh. + H.P. 
1 111 11 11 


-G i (i = 1, 2) 


(4.22) 


we obtain 

P T T T H.T. + T T H.T.P. = -tTg.T. (i = 1, 2) 
liii 1111 l l i 


(4.23) 


or by introducing 


H. = T T H.T. , G. = T T G.T. (i = 1, 2) 
l iii-i ill. ’ ’ 


we finally get 


P T H. + H.P. = -G. (i = 1, 2) 
ii ii l v J 


(4.24) 


(4.25) 


I i 

The transformations in (4.24) do not change the positive definiteness 

character of LL and CL and we arrive at Liapunov matrix equations (4.25) 

in transformed domain. Note, however, that in general, matrices r IL (i = 1, 2) 

T 

are not orthogonal since TLI\ f I and eigenvalues of G^ and G^ as well 
as of and FL (i = 1, 2) are distinct. Therefore, after solving the 

Problem C in transformed space, it is not possible to come back to the origi- 
nal description and obtain the same estimate of the coupling parameter. The 
whole procedure should be then carried out in transformed space. 

By retaining the cross -coupling terms only and neglecting self -coupling 
terms as the second order effect in (4.20) and (4.21) we arrive at: 
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VfLobtem C: 


± XCCh) X(H.) 

max y 

G A (H^ AfTT) 

subject to: P?lh + = ^ = ^ (4.26) 


Since the two problems (for i = 1, 2) are identical, we will how leave 
out index i for the simplicity of presentation and use P as a canonical 
representation for either subsystem. 

Before attempting to solve Vnoblrn C we will prove some preliminary re- 
sults . 

To be concise, we will define for the n x n matrix G = (g„) an 
p x p principal submatrix jas 


G p = G I 


1 1 1 2 X P V = 

i-i i 9 ■ • • i_. 

12 p. 


g i i g i i 

hi X 1 2 


1 2 1 1 1 2 1 2 


. 8 i i 

Pi P 2 


’i„i 


’Vp 


• g i 2 i 

2 p 


Vpj 


(4.27) 


where 1 < i. < i- . . . < i <n, and 1 < p < n . 

— 1 2 p — — — 

We will now state and prove the following theorem we need in order to 
solve Problem C and known as the Inclusion Principle [5] . 

Theorem 4.1 : Let X (G) and A(G) be the minimum and the. maximum eigenvalue. 

an n x n symmetric matfiix G and X (Gp) and A(Gp) the minimum and 
maximum eigenvalues any p x p principal AubmatAix Gp . Then ( {ok any 
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1 < p < n 


A (G) > A(G p ) 


X(G) < A (G p ) 


(4.28) 


pAQoj . First we recall the variational description of the minimum and maximum 
eigenvalue of a symmetric matrix G 

T 

, x Gx 

A (G) - max - rp - 


x^O x x 


T 

x x Gx 


A (G) — min — rjr 
x^O xx 


(4.29) 


where x T = (x^ , x^ , ... , x n ) . We obtain the vector x^ from the vector 

x by simple vanishing its components x^ = 0 for all i / ip i^, ... , i 

pT 

and 1 1 i 1 n , that is x r = (0, x. , . . . , x. , 0, . . . x. , . . . , 0) . 

X 1 1 2 x p 

Since the maximum on a set is no less than the maximum on any of its 
subsets, we have the following: 


‘ x T fx 

A(G) = max — — > max 


xpV 


xi<0 xV x p /0 x p ‘x P 


.T A ( G p) 


(4.30) 


Similarly the minimum on a set is no greater than the minimum on any of 
its subsets and therefore 


A (G) = min 
x^O 


T 

x*Gx 

T 


x x 


£ min 
xP^O 


xP T GxP 

T 

xP xP 


= A (G ) 
P 


(4.31) 


This completes the proof. 
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Note however, that this theorem is valid for more general case of Hermi- 

tian matrices, but for our purposes, it suffices to consider symmetric matrices 

since our whole analysis is in the real domain. 

As a consequence of the preceding theorem, we have the following: 

Let a symmetric matrix G be partitioned into blocks G = (G. .) , such 

k 

that G. . are n. * n. matrices, (i, j = 1, 2, . . . , k) , £ n- = n . De- 
n ii i=1 l 

note by G^ the block matrix obtained from G by taking G^ = 0 (zero 

matrix) for i f j . Then 


A(G) > max A (G i± ) = A^) , 1 < i < k (4.32) 

i 

A(G) <_ min A(G^) = A(G^) , 1 <_ i <_ k (4.33) 


Using this result we can state and prove the following theorem. 


Thz.on.rn 4.2 : When the. canonical nepn.es entation P ofi the subsystem n x n ma 

t/ux P with, complex eigenvalues , i = 1, 2, . . . , m and neal eigen 

values , i = m + 1 , . . . , k , such that crL < 0 , (i = 1 , ... , k) , n = 
m + k Is 


a l “l 
'“I G 1 


O 


P = 


a to 
m m 

-to a 
m m 


m+1 


O 


(4.34) 


then the solution G ofi the Vnoblem C is 
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o 



G = 2g 






(4.35) 


wheAe g aj> an aAbWiany fiexit poA'itive numbeA, and the. value o fa the maximum 
natio aj> 


1 

2 


\(G*) A (H*) _ 

9 

A (H ) 


min | a ■ j 
i 


1 < i < k 


(4.36) 


psioofi . In the proof we will again use "decomposition type" approach. First 
we split the matrix P into quasi -diagonal block form, that is, we represent 
P in (4.34) as P = (P„) , where 

for i = 1 , ... , m P^ = for i = m+1 , ... , k 
= 0 for i f j , (i, j = 1, 2, ... , k) (4.37) 



If we now use the same partitioning pattern for the matrices H and G 
from the Liapunov matrix equation 


P T H + HP = -G 


(4.38) 


we obtain the system of Liapunov-like matrix equations: 

^ii^ii + ^ii^ii = Ai 1 = 1 ’ 2 > *” ’ k 


(4.39) 
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P T .H. . + H. -P-- = -G.. i < j , (i, j = 1, 2, ... , k) (4.40) 

n 13 . 13 JJ ij 

The well known result that the matrix equations (4.40) have the unique 
solution for the unknown matrices is that A r (Pb ^) + A^ (Fb j ) f 0 for 

i < j , (i, j = 1, 2, ... , k) and r, i £ (1, 2} . This is really the case 
since the subsystem matrix P is supposed to be stable. Therefore, if 
Ghj = 0 , i < j , (i, j = 1, 2, ... ,k) then the unique solution of (4.40) 
is lb j = 0 , i < j , (i, j - 1,-2, ... , k) . 

Suppose now that we picked some G . We denote by G^ matrix obtained 
from G by letting CL. = 0 , i < j , (i, j = 1 , 2, ... ,k) , and similarly 
Hp obtained from H . Since e 0 implies ib j = 0 (for i / j) , then A 
using the results of Theorem 4.1, that is, (4.32) and (4.33) 

A(G) < \(G D ) 

A(H) < A(H D ) 

A(H)>A(H d ) (4.41) 


From (4,41) we obtain the following inequality for the ratio R(G) 


RCGj - ( 

‘ A(H) A (H) 


< 


1 A(G d ) A(H d ) 
I A(H d ) A(H d ) 


- 


(4.42) 


because G , H are positive definite matrices and therefore have real posi- 
tive eigenvalues. 

From (4.42) we conclude that the maximum of R(G) over G will be ob- 
tained for G Q , which implies . The proof is complete. 

The equations (4.40) are then automatically satisfied and we should solve 
the following subproblem: 
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Sub problem S: 


max; 

(G } 
P 


i HG p ) xoy 


I 71 


» (Hp) 


where P p H p * HpP p - -G p 


and 


P = 
P 


a a) 

P P 


-w a 
_ P P. 


Lama 4.3. The solution to the. AubpsiobZm S t& 


* 

G = g 
P P 


2 la | 0 

1 P 1 


0 2 la 


pU 


(4.43) 


(4.44) 


MkeAe. gp -u aAbttywA.y positive yumbeA and the maxA/num value ti 


VfLOoj. I£ we take 




the AubpAoblem S can be reformulated into 



/ 2 2 

The solution to (8.45) is obtained when y (g^- g^) + 4g 2 = 0 or when 
g^ = g^ and g^ = 0 , the maximum value is 
completed. 

Using the result of this Lemma we can now proceed to prove Theorem 4.2. 


and the proof of Lemma is 
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From (8-42) 


max R(G) = max R(G D ) = max y ^ 


1 -«h? 


D 


D 


Oh? 


, mm X(G..) mm X(H..) 
1 . v 11' • v n J 

l 


max y _i 

{G ii } [max A (H . . ) ] 2 


i=l, ,k 

The following inequalities can be derived from (4.46) 


li" 


(4.46) 


min X (G . . ) min X (H. . ) 

- if - 11 

max y - * < 

{G u ) [max A(H n )] Z 

i=l,2,...,k 


for every j = 1 , 2 , ... ,k. 



[ACH ^)] 2 


(4.47) 


Since the constraining equations (4.39) are completely decoupled and us- 
ing the result of Lemma 4.3 


max 



1 

2 


>( V 

tAffl ^)] 2 



j = 1, 


(4.48) 


However, inequality (4.47) is satisfied for every j = 1, 2, ... ,k, and 
therefore also for one which has minimum module of a real part of eigenvalue. 
Therefore , 

max R(G n ) < min ja.j 1 < j < k . (4.49) 

e D J ' 

If we now can find such a G^ that 

R(G D ) = min |a^| 1 1 j < k (4.50) 

we have found a solution to Problem C. 

This is, however, easy to accomplish simply by taking G^ = (G^) where 
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G. . 
11 


for all i = 1, 2, 


t 

to 

Q 

0 

~ * 

, which leads to H. . = 

’ XI > 

"1 

0 

0 

2 1 0 . | 

0 

1 






... , k . Therefore H is an n * n identity matrix with X (H) = A (H) =1 , 
XfGp) = min |a^| (j = 1, 2, . . . , k) , and equality (4.50) is satisfied. 

Note that G = g • Gp (g > 0) is also a solution to the Problem C 
since the ratio R(G) is invariant under multiplication of a chosen G with 
the positive scalar. This completes the proof of Theorem 4.2. 

But, still one question remains unanswered. That is, whether the solu- 
tion to Problem C in the transform domain gives the greater value for the maxi- 
mum allowable coupling parameters than the solution to the Problem C in the 
original domain? As it was noticed earlier, no direct correspondence has been 
established between the values R(G) and R(G) . Since G = T T GT , and T non- 
singular transformation matrix in general not orthogonal, T T f I , it was only 

certain that these values were different. One might naturally ask the question 

* 

whether there exists such a G in original domain which would give greater 

£ ft 

value for the ratio R(G ) than the maximum value of the ratio R(G ) = min|a- 1 

j 3 

(j = 1, 2, — , k) obtained in the canonical domain. We will answer this ques- 
tion negatively by proving the following theorem. 


Theofim 4.4: Let G , G 

, and H , H , 

be 

n x n Ayrmetsitc positive defitntte 

matAsLcet) AeZated by 




G = T T GT 

, H = T T HT 



T 

~T~ 



P H + HP = -G 

, P H + HP = 

-G 


P = T _1 PT 

, T T T f I 


(4.51) 


wheAe P -u the canonical siepAeA entatxo n the known Mxbt>yt>tm matsux P , 
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and T- -ci co k/lu ponding noAmallz&d Voavu^^onmation. Then 


R(G) = j X - ( % — < \ X(G) ? ^ = R(G) 
Z A Z (H) Z A Z (H) 


(4.52) 


Pfioo^ . Using the variational description o£ the minimum and maximum eigen- 
value we have 


x T Hx _ ic A T 

m — ■“ H13X f-jp rp . . 

~ ~ Cur> xj^O X T Tx 


A(H) = max 

x^O xx 


~ 1 'T 1 r r C C-Jn v T^T 


x^O x T Tx 


(4.53) 


since x = Tx and T is nonsingular. 

Now, from (4.53) we derive the following weak (greater than or equal) in- 
equality 


~T-~ 

X Hx 

A(H) = max -rp-K, — = max 
x^O • x J T Tx x^O 


~T~~ 

xHx 

-T- 
x x 

-T T - 
xTTc 

-T~ 
x x 


max*^ 

xj^O xx • 

- T - 
xT J Tx 


rain — ™ — 
#0 xx 


(4.54) 


T -t T- T 

Since x x is positive definite quadratic form so is x T Tx and T T is 

symmetric positive definite matrix. Furthermore, due to normalization of 

T = (t-^ t 2 ,...,t n ) where t^(i = 1, 2, ...» n) are column vectors, we have 

t^t- =1 i = 1, 2, • ,n . 

l l ’ ’ 


Then 


T 

T T = 


T 

+■ t 
J '2± 


T 

t t, 
n 1 


T T 
ti to tit 


12 13 

1 


T 

t t 

l 2 3 


T 

Vn 
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T 

Now, since T T f I at least some of the off-diagonal entries will be 
non- zero and direct application of the Theorem 4.1 gives 

A(T T T) < 1 

A(T T T) >1 (4.55) 

Using (4.55) we get 

A (H) > AQP > A(H) (4.56) 

~ l(TT) 

Similarly, for the minimum eigenvalue 

, T __ max 

A (H) = min ^ < A (H) (4.57) 

x^O x T Tx x T Tx A(T X T) 

ITlcLX Fp 

x/0 xx 

By the same procedure 

A (G) < A (G) (4.58) 

Now from (4.56), (4.57) and (4.58), the inequality (4.52) directly fol- 
lows and the theorem is proved. 

T 

Note, however, that in the very special case when T T = I then R(G) = 
R(G) . 

Stated in words, for every choice of G in original space we can find 
G in canonical space which would give greater (or in the very special case 
equal) value for ratio R(’) . Consequently, it must be true also for the 
maximum values, that is 

max R(G) < max R(G) (4.59) 

G G 
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* 

The inequality (4.59), however, does not imply that the matrices G and 

'J' 

G for which these maximums are obtained are related by G = T GT . There- 

* * 

fore, how to find G and R(G ) remains an open question, though, because 
of (4.59) we are not interested anymore in answering that question. 

As compared to (4.4) we have the comparison matrix A in the canonical 
domain 


A = 


~ x X(G X ) 

2 a(h 2 ) 


21 1 * 1 

^ (H 1 ) X(H 2 ) 



a(h 2 ) 

/ A(H 1 )A(H 2 ) 


1 *CG 2 ) 

2 A(H 2 ) 


(4.60) 


where and £ 2 1 are norms t ^ ie transformed intercoupling matrices Q 22 

and Q respectively. Due to the sparse structure of the coupling matrices 
Ql 2 ’ ^21 ancl taking normalized transformations T^ , T 2 we obtain 


^12 ? 12 » 5 21 < C 21 

what the computer results confirm. 

Similarly to (4.5), the second stability inequality gives 

2 ]_ X(G 2 ) A(H 2 ) A(G z ) A(H 2 ) 

J i 12 i 21 a 2 (h 2 ) A 2 (H 2 ) 


(4.61) 


(4.62) 


If we now choose and C 2 as proposed by Theorem 4.2 and use the re- 
sult of Theorem 4.4, that is (4.59), together with (4.61), we find that the 
maximum stability region of the coupling parameter £ is obtained for the 
comparison matrix 
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An 

_5 2i Ul " 1 0 m2 ^ 

where o and o are the eigenvalues with the minimum module of the real 
part for the subsystem matrices and . The largest estimate for the 
stability region 

(4.64) 

is obtained for the subsystem Liapunov functions 

V 1 = I 1*1 1 I ’ v 2 = 1 Ix 2 1 I (4.65) 

Before presenting the numerical results which will show to what extent the re- 
sults have been improved by this optimization approach, we will discuss some 
of the features introduced by the canonical transformations. 

The crucial advantages of introducing the canonical transformation lie in 
the following: 

1. The problem which did not seem to be tractable appeared in a 
manageable form in the canonical space, where it was 'Veil 
defined" and solved. 

2. Due to the relationships (4.59) and (4.61), the solution in 
the canonical space offered the solution to the original prob- 
lem of maximizing the stability region of the parameter ? . 

3. The canonical transformation is introduced at the subsystem 
level. The subsystems are assumed to be of the low order. 
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The numerical problems imposed in canonical transformation 
are then of much less difficulty than numerically solving 
the problem using nonlinear mathematical programming algor- 
ithms. 

4. The simple analytical expression (4.64) not only offers the 
practical way of computing maximum allowable 5 but more 
importantly, gives the significant insight into the exist- 
ing trade-off between the degree of the subsystem stability 
and the strength of interconnections as compared to the over- 
all system stability. This immediately suggests a possible 
design procedure based on introducing the controls on the 
subsystem level in such a way as to shift all subsystem 
poles as much to the left but without too much affecting 
the interconnections. These qualitative discussions and 
the existing trade-off will be made apparent in the numeri- 
cal examples of Passive and Active Control case. 

The computer program implementing this optimization procedure is given in 
the Appendix together with the obtained results for both Passive and Active 
Control case. The results are as follows: 


PASSIVE CONTROL CASE 


The two subsystem matrices are 


P l = 


0 

0.8478 

0 

0 


0.0463 

0 

0 

-1.0463 


0.4387 

0 

0 

-3.2253 


0.0131 

0 

1.0 

-0.0665 
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0 

1.0 

0 

0 


-1.7865 

-0.0534 

0 

2.0 

P 2 = 

0 

0 

0 

1.0 


0 

-2.0 

-904.5314 

-1.2036 


The canonical forms are 


p i* 


-0.0005 

-0.2909 

0 

0 


0.2909 

-0.0005 

0 

0 


0 

0 

-0.0327 

-1.7607 


0 

0 

1.7607 

-0.0327 


P 2 ' 


-0.0266 

-1.3334 

0 

0 


1.3334 

-0.0266 

0 

0 


0 

0 

-0.6019 

-30.1359 


0 

0 

30.1359 

-0.6019 


The computed norms of the transformed intercoupling matrices are 
l u = 0.6314 < 0.7766 = l u 
l 21 = 1.4545 < 1.8478 = l n 

The comparison matrix is 


-0.0005 0.6314jc| 

A = 
m 

1.4545 | 5 | -0.0266 

and the largest estimate of the stability region of £ 


(4.66) 


(4.67) 


(4.68) 


(4.69) 
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(4.70) 


|5| < 5 = 0.403 x 10 


ACTIVE CONTROL CASE 


The subsystem matrices are 


h - 


0 

0 

0 

0.0463 

0 

1.0463 


0 

0 

0 

-1.2456 

-0.8478 

-1.2456 


0 

0 

0 

-0.4387 

0 

-3.2253 


1 

0 

0 

CL. 2456 
-0.1521 
-1.2456 


0 

1 

0 

1.0463 

0 

1.0463 


0 

0 

1 

-0.0131 

0 

-0.0665 


P 2 


-0.1823 

0 

0.1823 

0 

-2 


0.0642 

0 

-1.8508 

0 

0 


0.0019 

1 

-0.0553 

0 

-2 


0 

0 

0 

0 

-904.5314 


0 

0 

2 

1 

-1.2036 


The canonical forms are 


h- 


-0.1288 

0 

0 

0 

0 

0 


0 

-0.8408 

0 

0 , 
0 
0 


0 

0 

-0.0783 

-1.7024 

0 

0 


0 

0 

1.7204 

-0.0783 

0 

0 


0 

0 

0 

0 

-0.0929 

-0.8990 


0 

0 

0 

0 

0.8990 

-0.0929 


(4.71) 
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P 2 = 


-0.1761 

0 

0 

0 

0 


0 

-0.0306 

-1.3567 

0 

0 


0 

1.3567 

-0.0306 

0 

0 


0 

0 

0 

-0.6019 

-30.1359 


0 

0 

0 

30.1359 

-0.6019 


(4.72) 


The computed norms of the transformed intercoupling matrices are 


? 12 = 14.2154 < 314.5269 = l u 


(4.73) 


but 


l 2 i = 3.7497 > 2.7088 = £ 21 


? 12 ‘ C 21 < 5 12 5 21 


(4.74) 


The comparison matrix is 


A = 
m 


-0.0783 


3.7497UI 


14.2154 5 I 


-0.0306 


(4.75) 


and the largest estimate of the stability region of £ is 

|?| 1 S m =. 0.6714 x 10" 2 (4.76) 

If we now compare the results of direct approach taken in reference [1] 

\ 

with the optimization approach, in both the Passive and the Active Control case, 
we find considerable improvement of more than 10^ times. Besides that, random 
experimentation with different choices G^, G 2 in the case without optimization, 
showed that chosen identity matrices for G^ and G 2 appear to be the best. 
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These results reveal the significance of introducing the canonical trans- 
formation on the subsystem level. Besides the significantly improved results, 
the computational alternative offered by this optimization procedure seems to 
be very attractive and rather simple. There is no need for solving subsystem 
Liapunov matrix equations and computing the eigenvalues of the G , H matrices. 
These computational steps are replaced by the subsystem canonical transforma- 
tions, whereas the other steps of computing the norms and the estimate of the 
stability region of the coupling parameter £ remain unchanged. 

The remarkable advantage of this approach is that the optimization problem 

does not require the numerical optimization procedure over all possible , 

G^ , but is based on the similarity transformation of subsystem matrices to the 

properly chosen canonical form. Therefore, the low-order of the subsystems, 

the simplicity of the optimization procedure, and substantial enlargement of 

4 

the stability region (more than 10 times) offer very attractive ways to cir- 
cumvent the inherent over-sufficiency of decomposition-aggregation stability 
procedure via the Liapunov vector functions. 

Furthermore, in the following section, these theoretical as well as compu- 
tational results will be directly applied to obtain optimal comparison system 
as suggested at the beginning of this section. 
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5 . OPTIMUM COMPARISON SYSTEM 


As we were motivated at the beginning of the previous section, we are 
really looking for the comparison system as given by (5.2). Before trying 
to formulate the optimality conditions for obtaining optimal comparison sys- 
tem, we will investigate some of its peculiarities . 

Suppose that we have the following comparison system 


V 1 1 'Vl + a l2 V 2 


V 2 — a 21 V l ‘ a 22 V 2 


A = 


" a ll °12 


a 21 " a 22 


(5.1) 


where a_ (i, j = 1, 2) are nonnegative numbers and > 0 , v 2 _> 0 . Due 
to linearity of the comparison system (5.1) after introducing the new variables 


V 1 ’ v 2 b >’ 


V 1 " «l v l v 2 = «2*2 (E 1> *2 * 0) 

we obtain the comparison system 


(5.2) 


V 1 - - a ll v l + “12^2 


v 9 < — a 01 V. - a 00 V„ 

2 — g 21 1 22 2 


(5.3) 


g 2 

where g = — > 0 

The comparison systems (5.1) and (5.3) are equivalent in the sense that 
one may be obtained from the other and vice versa by using (5.2). In terms 
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of subsystem Liapunov functions it amounts to arbitrarily choosing scalar mul- 
tiples of some known v^ , v 2 as the new Liapunov functions. 

If we consider a ^j (i, j = 1, 2) as the adjustable parameters depending 
on the matrices as given in (-5.4), we can formulate the following 

criterion for the optimality of the comparison system (5.4). 

PAoblm V . 

Find A 0 such that 

A° - A < 0 , VA 6 K (5.4) 

subject to constraints: 


P T H. + H.P. = G. 
IX xi x 


i = 1, 2 


A(G X ) = 2 



X(G 2 ) - 2|c m2 | 


(5.5) 


where K is the set of all matrices {G^, G 2 > which are normalized by the 
condition (5.5). Here , a m2 are as defined by (5.63). 

Having in mind the equivalence between the systems (5.1) and (5.3) we 
can now reformulate this problem into 


VKoblm E. 


and 


max 


G. 

x 


1 

2 AIH7J 


i = 1, 2 


min 

{G 1’ G 2 } 


a ay a(h 2 ) 

ray a(h 2 ) 


subject to the constraints: 
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PTH. + H.P. = -G. i = 1, 2 
1111 l 


*( G l) 2 1 °mi i x( ^ G 2^ = 


(5.6) 


Following exactly the same steps as in solving Psioblm C and then Problem 
_ _ * ~ * 

C we arrive at the same solution , G^ for the ?no blm E. Therefore, the 
obtained comparison matrix A in (4.63) is also the optimal comparison matrix 
A° . 

Furthermore, with this kind of approach there is no more need to neglect 

the self-coupling terms as was done in arriving at (4.4) but the complete sys- 

£ & 

tern (4.20), (4.21) can be optimized by the same optimal choice of G^ , G^ . 
This is true, because additional optimality requirements 


mm 

G. 


A(H.) 

ITiTT 


i = 1, 2 


(5.7) 




are satisfied for the same optimal choice of G^ , G 2 . 

Consequently, the optimal comparison matrix for the Passive Control case, 
taking into account the self -coupling terms, is 


A° = 


-0.0005 + 0.7665 


1.454k 


0.631k 


-0.0266 + 0.144?' 


(5.8) 


and for the Active Control case is 


A° = 


-0.078 + 0.1655 


3.749 5 


14.215 1 5 I 


-0.030. + 0.7445 


(5.9) 


Finally, we have succeeded in obtaining the optimal comparison systems for 
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the mathematical model description of the Passive and Active Control Case of 
the Spinning Skylab and at the same time we have maximized the stability re- 
gion of the important structural parameter 5 which represents the ratio of 
the asymmetry of the boom settings and the length of the booms. 
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6. CONCLUSION 


A decomposition-aggregation method is outlined for stability analysis 
of large-scale dynamic systems. The method takes advantage of special struc- 
tural features of the complex systems to reduce the memory and computational 
time requirements. Furthermore, the method is suitable for accommodation of 
nonlinearities either in the subsystems or in their interconnections and pro- 
vides the significant insight into the structural properties of complex dyna- 
mic systems . 

The straightforward application of the Liapunov Direct Method to the 
stability analysis of high order multivariable dynamic systems though theo- 
retically possible, involves formidable practical limitations. By decompo- 
sition of the system into a number of lower order subsystems, it has only 
been necessary to apply Liapunov Direct Method to the lower order subsystems. 
Therefore, the full advantage of the experience already accumulated in the 
construction of Liapunov functions for systems of comparatively small order 
can be taken. Consequently, this approach broadens the class of problems to 
which Liapunov Direct Method has practical application and also makes the 
treatment of large-scale systems feasible. However, the method is inherent- 
ly conservative since a series of approximations are involved in obtaining 
the comparison inequalities and establishing sufficient conditions for sta- 
bility. This may lead to overly restrictive (overly sufficient) conditions' 
for stability as the results in reference [1] point out. 

The decomposition-aggregation method was applied in reference [1] to the 
dynamic model of a spinning Skylab. After the model was decomposed into the 
wobble and spin subsystems, both the passive and the active control were con- 
sidered. Such decomposition made an important structural parameter asymmetry 
in the boom settings to appear as an interconnection parameter of the two sub- 
systems. Subsequent stability analysis was aimed at estimating the interval 


50 



of this coupling parameter for global exponential stability of the overall 
system. The estimates obtained in reference [1] turned out to be very con- 
servative since the flexibility of the decomposition-aggregation method was 
not used to the full extent. 

In this report, the problem of maximizing the stability region of the 
coupling parameter was formulated and resolved as a well-defined optimiza- 
tion problem [4] . The obtained optimality conditions turned out to be ra- 
ther simple, requiring only the canonical transformation of the subsystem 
matrices as well as corresponding transformation of coupling matrices. More- 
over, the possibility of introducing the controls at the subsystem level and 
still maintaining the insight into the effect on the overall stability region, 

provides a basis for possible design procedure. The substantial enlargement 

4 

of the stability region (more than 10 times compared to nonoptimal procedure 
outlined in reference [1]) and the simplicity of the optimization procedure 
offers a very attractive way to overcome the inherent over -restrictiveness of 
the decomposition-aggregation stability procedure via the Liapunov vector func- 
tions. As shown in the Appendix, the entire optimization procedure can be 
suitably programmed for computer applications. 
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APPENDIX 
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DESCRIPTION OF THE COMPUTER PROGRAM 


The program for maximizing stability region of the coupling parameter 5 
and obtaining the optimum comparison system is realized through the canonical 
transformation of the subsystem matrices to properly chosen canonical form. 

The method used to compute the canonical transformation matrix T and 
reduce matrix P to its canonical form P as given by (4.34) is as follows: 
The coefficients of the characteristic equation are computed using 
Leverrier algorithm. Then, the roots of the characteristic equation are found 
using standard subroutines and matrix P is formed. The system of linear 
equations 

PT - TP = 0 (A.l) 

is then transformed into 

[P® I - I® P] t = 0 (A. 2) 

where the unknown vector t is formed from the columns of the matrix T and 

® denotes Kronecker product. Since N x N matrix T is not unique, first 

2 2 

row is chosen arbitrarily and the corresponding system of (N - N) * (N - N) 
linear equations is solved to obtain the transformation matrix T . 

Besides authors' own subroutines whose descriptions are presented here, 
two other IBM supplied subroutines were used: SIMQ for solving the system of 

linear equations and POLRT for finding the roots of the polynomial. 
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Subroutine. ENORM 


Purpose: Compute the transformed coupling matrix and its norm. 

linage: CALL ENORM (Q, T, TR, K, L, TINV, SIM, TEMP, EST) 



Q - Original coupling matrix 

T - First canonical transformation matrix of order K x K 
TR - Second canonical transformation matrix of order L x L 
K - Number of rows in matrix Q 

L - Number of columns in matrix Q 

TINV - Inverse matrix of the matrix T 

-1 T -i 

SIM - Vector of eigenvalues of the matrix (T QTR) (T QTR) 

TEMP - Temporary storage matrix 

EST - Norm of the transformed coupling matrix 

Remark : The computed transformed coupling matrix T ^QTR is stored in Q . 

Method : The norm of the matrix Q is- computed as the square root of the 

maximum eigenvalue of the matrix QTq . 

Subroutine CHVAL 

Purpose: Compute the characteristic polynomial and the characteristic 

values of the given matrix. Then form the canonical matrix as 
given by (4.34). 

th> qq e •' CALL CHVAL (A, N, XCOF, COF, ROOTR, ROOTI, B, C, IER) 

Ve&cJviption ofa parameters : 1 

A - Name of input matrix of order N x N ■ 

•A 

N - Order of input matrix 

XCOF - Vector of N + 1 coefficients of the characteristic polynomial 
of matrix A ordered from smallest to largest power 

COF - Working vector of length N + 1 

ROOTR - Vector containing real parts of eigenvalues of A 

ROOTI - Vector containing imaginary parts of eigenvalues of A 



B - Canonical matrix of a matrix A 

C - Temporary storage matrix 

IER - Error code where: 

IER = 0 No error 

IER = 1 N less than one 

IER = 2 N greater than 36 

IER = 3 Unable to determine root with 500 
iterations in 5 starting values 

IER = 4 High order coefficient is zero 
SabAoivtinz CRONE 

PuApo-i e •• From the matrices P and PS form the matrix P©I - I ©PS where 

@ denotes Kronecker product . 

Uiages CALL CRONE (P, PS, CRON, N) 

VoAcfyiption o{, paArnieXeAA : 

P - Name of one input matrix 

PS - Name of second input matrix 

CRON - Name of output matrix P®I - I ©PS of order N 2 x N 2 
N - Order of matrices P and PS 

Subfiouutim SLINE 

PuApo4e.* Find canonical transformation matrix by solving system of linear 
equations PT - TPS = 0 where P and PS are known N x N 
matrices . 

U6age: CALL SLINE (CRON, XR, XC, N, KS) 

* 

Via oviption oj paAamiiteAA : 

CRON - Input matrix of the dependent system of homogeneous linear equa- 
tions PT - TPS 

XR - Input vector consisting of elements of the first row of transfor- 
mation matrix T 

XC - Computed vector of coefficients of the independent system of 

(N 2 - N) x (N 2 - N) linear equations obtained from PT - TPS = 0 
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N 


Order of matrices P and PS 


KS - Output digit 

0 for a normal solution 

1 for a singular set of equations 

2 2 

Method: Divide N * N matrix CRONE into 4 block matrices where diagonal 

2 • 2 

block matrices are of order N x N and (N - N) ’ (N - N) . Use 

then second one as the coefficient matrix of the system of linear 

2 

equations. Multiply then (N - N) x N off-diagonal block matrix 

with the known N * 1 vector XR to compute the coefficient vec- 
2 2 

tor of the (N - N) * (N - N) system of linear equations. 
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PAGE 

1 

09/25/73 

JOB 

NO. 0AAI03 

// JOB 

T 

00 Ml 



LOG DRIVE 

CART SPEC 

CART AVAIL 

PHY DRIVE 

0000 


0001 

0001 

0000 

V2 M 10 


ACTUAL 16K 

CONFIG 1 6K 



// ♦ 

// ♦ PASSIVE AND ACTIVE STABILITY OF THE SPINNING SKYLAB 

// * M A X IMA 2 1 NG THE REGION OF STABILITY WITH RESPECT TO COUPLING PARAMETA, 

// ♦ 

// FOR 

♦ONE WORD INTEGERS 
♦CIST SOURCF PROGRAM 

SUBROU T i NE ENORM < Q , T , TR , K » L , T I NV,.S I M , TEMP., E ST.) ... . 

DIMENSION LL ( 12 ) ,MM{ 12) 

DIMENSION Q ( 1 ) , T ( 1 ) , TR ( 1 > , T I NV< 1 ) ,SIM(1) ,TEMP(1) 

CAU. MCPYIT »TINV»K,K,0) 

CALL MINV(TINV»K,OET»LL»MM) 

CALL G v iPRD(TINV,Q,TEMP,K,K,L) 

CALL C LIP R D ( T E M P , T R , Q , K , L , L_) 

CALI. MATA(Q,SIM,L,L»0) 

CALI. F 1 GEN ( S IM, R , L ,1 ) 

EST---SDRT(SIMm) 

RETURN . 

END 


FEATURES SUPPORTED 
ONE WORD INTEGERS 

CORE REQUIREMENTS FOR ENORM 

COMMON 0 VARIABLES 30 PROGRAM 98 

RELATIVE LNTRY POINT ADDRESS IS 0020 (HEX) 

END OF COMPILATION 

II DUP 


♦STORE WS UA ENORM 

CART ID 0001 DB ADDR 5BB0 . DB_ CNT 0007 

// EJECT .... ; ..... . 
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PAGE 


2 


09/23/73 


// FOR 

♦ LIST SOURCE PROGRAM 
♦ONE WORD INTEGERS 

SUBROUTINE CHVAL < A , N , XCOF ,COF , ROOTR , ROOT I , B , C , I ER) 
DIMENSION All), B(l), Cil) ,XCOF(l) , C0F(1) ,R00TR(1) , ROOT 1(1) 
I SUB ( I , J ) = ( J-1)*N+I 

\m=n*n 

DO 10 I = 1 » NN 

10 B( I 1=0.0 
DO 12 1=1, N 
. J=I+( 1-1 )*N 
12 B ( J ) = 1 . 0 

DO 20 I = 1 , N 

CALL GMPRDI A,8,C,N,N,N) 

• TRACE=0. 0 
DO 15 K= 1 , N 
J=K+(K-1)*N 
15 TRACE=TRACE+C(J) 

ROOTR ( I ) =TRACE 

20 CALL MCPY(C,B,N,N,01 
COF( 1) = 1.0 
DO AO K = 1 , N 
SUM=0. 0 

00 30 1 = 1, K 

_ L = K- I + 1 

30 SUM=SUM+COF(L)*ROOTR< I ) 

AO COF ( K+l 1 =-SUM/K 

M = N+ 1 

DO A 5 1=1, M 
J =N- I +2 

A 5 XCOF ( I )=COF( J) 

CALL POLRTt XCOF, COF, N, ROOTR, ROOT I , IER) 

DO 50 1=1, N 
50 ROOT I { I 1 =-ROOT I ( I ) 

DO 60 1=1, N 
DO 60 J = 1 , N 

. JJ=ISUB( I , J ) 

60 B ( J J 1=0.0 
DO 70 1=1, N 

1 I = I SUB ( 1,1) 

D ( I I 1 =ROOTR ( I 1 
I F ( ROOT 1(1)1 71,70,73 

73 I J = I SUB (1,1 + 1) 

B ( I J 1 =ROOT I ( I 1 
GO TO 70 

71 I J = I SUB (1,1-1) 

B ( I J 1 =ROQT I ( I 1 
70 CONTINUE 
RETURN 
END 

FEATURES SUPPORTED 
ONE WORD INTEGERS 

CORE REQUIREMENTS FOR CHVAL 

COMMON 0 VARIABLES 20 PROGRAM 528 
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PAGE 


3 


09 / 23/73 


RELATIVE ENTRY POINT ADDRESS IS 002.D (HEX) 

END OF COMPILATION 
// DUP 

*STORE WS UA CHVAL 

CART ID 0001 DB ADDR 5BB7 DB CNT 0023 

// EJECT ■ 
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PAGE 


4 


09/23/73 


// FOR 

♦ LIST SOURCE PROGRAM 

♦ ONE WORD INTEGERS 

SUBROUTINE CRONE ( P , PS , CRON , N ) 

DIMENSION PtI),PS(l) » CRON 1 11 

I SUB ( I , J )=( J-l )*N+I 

ISUBC(I.J)=(J-1 )*N*N + I 

DO 20 K= 1 » N 

DO 20 L = 1 » N 

I P= I SUB ( K » L ) 

IF(K-l) 22,21,22 . 

21 DO 23 1 = 1, N 

DO 2 3 J = 1 » N 

I I = { K- 1 ) *N + I 
' JJ=(L-1)*N+J 

NN= I SUBC ( I I , JJ) 

IPS=ISUB(J, I ) 

IF(I-J) 24,25,24 
... 25 CRON ( NN ) =P ( IPl-P SI I PS) 

GO TO 23 

24 CRON(NN)=-PS( IPS) 

23 CONTINUE 
GO TO 20 

22 DO 26 1=1, N 

DO 26 J = 1 , N 

I I = ( K- 1 ) #N+ I 
JJ»( L-1>*N+J 
NN = l SUBC ( I I , JJ) 

IF(I-J) 28,27,28 ... 

27 CRON ( NN ) =P ( IP) 

GO TO 26 

28 CRON ( NN ) =0. 0 
26 CONTINUE 

20 CONTINUE 
RETURN 
END 

FEATURES SUPPORTED 
ONE WORD INTEGERS 

CORE REQUIREMENTS' FOR CRONE _ 

COMMON 0 VARIABLES 12 PROGRAM 312 


RELATIVE ENTRY POINT ADDRESS IS 0036 (HEX) 


END OF COMPILATION 
// DUP 

♦STORE WS UA CRONE 

CART ID 0001 . OB AODR 5BDA OB CNT 0014 
// EJECT 


61 
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// FOR 

♦ ONE WORD INTEGERS 
♦LIST SOURCE PROGRAM 

SUBROUTINE SL I NE ( CRON , XR , XC , N , KS ) 
DIMENSION CRON ( 1 ) » XR ( 1 ) » XC ( 1 ) 

I SUB ( I , J ) = ( J-n*N+I 

NN=N*N __ 

M = N* ( N-l ) 

MN=M*N 

KK*MN«(N-1) 

DO 10 J = 1 , NN 
DO 10 1=1, M 

N I = ( J- 1 ) ♦NN+N+I 

NJ=( J-l >*M+I 
10 CRONINJ)=CRON(NI ) 

CALL GMPRD(CR0N,XR,XC,M,N,1) 

DO 20 1=1, M 
20 XC I I ) =-XC( I ) 

00 30 1=1, KK 

L L = M N + 1 

30 CRON ( I I =CRON ( L L ) 

CALL SIMQ(CRON,XC,M,KS) 

DO AO 1 = 1, M 
K=NN- I + 1 

L = M- I + 1 

AO XC ( K ) =XC { L ) 

DO 50 1=1, N 
50 XC ( I ) =XR ( I ) 

DO 60 1=1, N 
DO 60 J = 1 , N 

1 F ( I -J ) 61,60,60 

61 I I = I SUB ( I, J) 

j j = i sub ( j , n 

•TEMP = XC ( II ) 

XC ( I I >=XC( JJ ) , 

XC(JJ)=TEMP 

60 CONTINUE 

RETURN 

END . . 

FEATURES SUPPORTED 
ONE WORD INTEGERS 


CORE REQUIREMENTS FOR SLINE 
COMMON 0 VARIABLES .. 18 PROGRAM . .368 

RELATIVE ENTRY POINT ADDRESS IS . 0025 (HEX) 


END OF COMPILATION 


// DUP . 

♦STORE WS UA SLINE 

CART ID 0001 DB ADDR 5BEE DB CNT 0018 


// EJECT 


62 
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// FOR 

*L I ST SOURCE PROGRAM 
* I OC S ( C ARD » 1403 PRINTER 
*ONE WORD INTEGERS 

REAL II , 12, I3,MASS,K1,K2 

DIMENSION P( 36) ,PS( 36) ,CRON( 1296) , XR ( 6 ) , XC ( 36 ) , TEMP ( 36 ) , T I NV { 36 ) 
DIMENSION T ( 36) ,Q{ 36 ) , SIM( 2 1 ) ,COF ( 7J , XCOF ( 7 ) , ROO TR ( 6 ) , ROOT I ( 6 )._ 
DIMENSION SM IN ( 2 ) 

C READ THE PHYSICAL PARAMETARS OF THE SKYLAB 

RE AD (2, 50) I 1, 1 2,13, G2, MASS, EK1 ,EK2, OMEGA 
50 FORMAT ( 8F 1 0 . 0 ) 

C COMPUTE THE NORMALIZED PARAMETARS 

K 1 = ( 12-13) /I 1 
K2= { I 3- I 1 ) / I 2 
ALPHA=( 1 . 0 + K 1 ) / ( 1.0-K2) 

BETA* I 3/ 1 1 

GAMA=(2.*MASS*G2*G2 ) /ll 
SIGSl=EKl/( MASS* QMEGA*OM EGA) 

SIGS2 = EK2/( M ASS*OMEGA*OMEGAJ_ 

SIGM1=SQRT(SIGS1) 

SIGM2=SQRT( SIGS2) 

DEL1=0.04*SIGM1 
DEL2=0.04*SIGM2 
GAMA 1 = 1. O-GAMA 

...... GAMA2=GAMA1*GAMA1. 

GAMA3=GAMA/BETA 

C 

C SUBSYSTEM ANALYSIS 

C 

DO 83 K= 1 , 2 

C READ N - THE ORDER O F T HE SUBSYSTEM 

READ (2,100) N 

100 FORMAT ( 12) 

NN = N*N 
NNN=NN*NN 

101 FORMAT ( 8F 10 . 0 ) 

DO 86 J = 1 , N N 

86 P(J)=0.0 

GO TO (41,42) , K 

C COMPUTE THE ELEMENTS OF THE FIRST SUBSYSTEM MATRIX PI 

41 P ( 5 ) = ( K 1 +GAMA ) / GAMA 1 

P ( 9) = ( GAMA*S IGS1 ) /GAMA1 

P( 13)=(GAMA*DELI)/GAMA1 

P ( 2 ) =K2 
P ( 1 5 ) = 1 . 0 

P ( 8 ) =-( Kl+1 .0) /GAMA1 
P( 12 ) = (GAMA-SIGSl-l.O) /GAMA1 
P( 16)=-DEL1/GAMA1 
GO TO 8 

C COMPUTE THE ELEMENTS OF THE SECOND SUBSYSTEM MATRIX P2 

42 P{ 5) =1.0 
P(2)=-SIGS1 
P ( 6 ) =-DEL 1 
P( 14)=2.0 

P ( 1 5 ) = 1 . 0 

P(8)=-2.0 

P( 12)=1.0-SIGS2 

P ( 1 6 ) =-DEL2 co 
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C W RITE The subsystem matrix 

8 WRITEI5, 160) 

160 FORMAT ( ' 1THE SYSTEM MATRIX IS*,/) 

DO 6 1 = 1, M 

6 WR I TE ( 5 , 200 ) < P ( J ) ,J=I , NN » N ) 

200 FORMAT ( IX, 12F 10. 4 ) 

C COMPUTE THE CHARACTERISTIC EQUTION AMD ITS ROOTS . 

CALL CHVAL(P»N,XCQF*COF, ROOTR , ROOT I , PS , TEMP , I ER ) 
WRITE!5,105) 

I F ( IER-3) 60,49,60 
49 WR I T E ( 5 , 300 ) 

300 FORMAT ( ' UNABLE TO DETERMINE ROOT WITH 500 ITERATIONS') 

GO TO 90 

60 WRITE! 5,400) 

400 . FORMAT ( • OTHE ROOTS OF THE CHARACTERISTIC EQUATION ARE',/) 
WR I T E ( 5 , 450 ) 

450 FORMAT!' REAL PART », 5X ,' I MAGI NARY PART ',/ ) 

DO 80 1=1 ,N 

. . 80 WR I T E ( 5 , 500 ) ROOTR ! I ) , ROOT I. (.1 .)_ 

500 F0RMAT(F12.4,5X,F12.4,/) 

WR I TE ( 5 , 250 ) 

C WRITE THE SIMILAR MATRIX 

250 FORMAT!' THE SIMILAR MATRIX IS',/) 

00 40 1=1, N 

...40 WRITE (5, 200) ( P S ( J ) , J = I , NN » N ) _ 

C FIND THE ROOT WITH THE MAXIMUM REAL PART 

SMIN(K)=R00TR(1) 

DO 81 1=1, N 

IF(SMIN(K)-ROOTR( I ) ) 82,81,81 ' 

82 SM IN ( K ) =ROOTR { I ) 

81 CONTINUE 

C READ THE VECTOR XR TO DETERMINE THE UNIQUE 

C SIMILARITY TRANSFORMATION 

RE AD ( 2 , 1 0 1 ) XR 
CALL CRONE (P, PS, CRON, N) 

WRITE(5, 105) 

. -.105 FORMAT!///) 

C SOLVE THE SYSTEM OF LINEAR EQUATIONS AND 

C COMPUTE THE SIMILARITY TRANSFORMATION MATRIX 

CALL SLINE(CRON,XR,XC»N,KS) 

I F(KS-l) 32,31,32 

31 WRITE(5,201) 

201 FORMAT!* SINGULAR CASE'..) 

GO TO 90 

C WRITE THE SIMILARITY TRANSFORMATION MATRIX 

32 WR I TE ( 5 , 305 ) 

305 FORMAT! • THE SIMILARITY TRANSFORMATION MATRIX IS',/) 

DO 12 1 = 1, N 

12 WRITE (5, 103) ( XC.(_J )j J= I,NN ,_N ) ; 

103 FORMAT (10F12.4) 

WRITE(5,105) 

IF(K-l) 83,84,83 

C CHANGE THE ORDER OF THE. FIRST SUBSYSTEM TO VARIABLE M 

84 M = N 

DO 85 J = 1 , NN ._ 

85 T ( J ) =XC ! J ) 

83 CONTINUE 

64 . 


C 
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C__ I N T e R C O N . N E C ,T_I 0_N A NA LYSIS 

C 

MN=M*N . . 

MM=M*M 

DO 61 1=1, MM 

C COMPUTE THE TRANSFORMED SELFCOUPLING MATRIX QII 

61 Q ( I ) =0 . 0 

Q( 5) =GAMA* ( 2.0*GAMA+K1-1 .0) /GAMA2 
Q(9)=SIGS1*GAMA#GAMA/GAMA2 
Q( 13)=(DEL1*GAMA*GAMA)/GAMA2 
Q( 2 ) =GAMA*( K2+1.0) /ALPHA 
0(8)=— Q( 5) 

Q ( 1 2 ) =-Q ( 9 > . _ 

Q ( 16 ) =-Q (13) 

CALL ENORM(Q,T,T,M,M,.TINV, SIM, TEMP, EST1) 

WR I TE ( 5 » 30 1 ) 

301 FORMAT! • 1THE TRANSFORMED SELFCOUPLING MATRIX Qll IS',/) 

DO 21 1 = 1, M 

; 21 WRITE(5, 103) ( Q ( J ) ., J- I ,MM , MJ. 

DO 62 1=1, MN 

C COMPUTE THE TRANSFORMED. INTERCOUPLING MATRIX Q12 

62 Q! I ) =0.0 
Q(5)=2.0*GAMA/GAMA1 
Q(9)=-GAMA/GAMA1 

Q ( 2 ) =GAM A* ( SIGSl+1.0) /ALPHA 

Q( 6) =GAMA*DEL1/ALPHA 
0 ( 8 ) =-Q ( 5 ) 

Q ( 1 2 ) =-Q ( 9 ) 

CALL ENORM ( Q , T , XC , M , N , T I N V , S I M , TEMP , E ST2 ) 

WRITE! 5, 302) 

...302 FORMAT!///, '. THE. TRANSFORMED INTERCOUPLING MATRIX Q12 IS' ,/ > 
DO 23 1=1, M 

23 WRITE ( 5, 103) ( Q ( J ) , J = I , MN , N ) 

C COMPUTE THE TRANSFORMED INTERCOUPLING MATRIX 021 

DO 63 1=1, MN 

63 Q( I )=0.0 

Q ( 2 ) =- ! 1.0 + K2) 

Q(8)=(2.0*GAMA+K1-1.0)/GAMA1 
0! 12)=GAMA*SIGS1/GAMA1 
Q ( 16)=GAMA*DEL1/GAMA1 

CALL ENORM { Q , XC , T , N ,M, T INV, SIM , TEMP , EST3 ) . 

WRITE(5»303) 

303 FORMAT!///, • THE TRANSFORMED INTERCOUPLING MATRIX Q21_IS'.,/) 
DO 25 1=1, N 

25 WRITE! 5, 103) ( Q ( J ) , J= I , MN , M ) 

C COMPUTE THE TRANSFORMED SELFCOUPLIN MATRIX Q22 

DO 64 I=1,NN 

64 Q( I ) =0.0 

Q( 2) =-GAMA*( SIGS 1+1.0) /ALPHA 
0(6) =-GAMA*DELl /ALPHA 
Q(8)=2.0*GAMA/GAMA1 ' 

Q ( 12 ) =-GAMA/GAMAl 

CALL ENORM ( Q, XC, XC, N, N,T INV, SIM, TEMP, EST4) 

WR I TE ( 5 , 304 ) 

304 FORMAT!///, • THE TRANSFORMED SELFCOUPLING MATRIX Q22 IS',/) 
DO 27 1=1, N 

27 WRITE(5,103) ( Q ( J ) , J= I , NN , N ) 

WR I TE ( 5 , 306 ) 
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306 FORMAT (////»' THE FOUR NORMS OF THE COUPLING MATRICES ARE',/) 
C WRITE THE FOUR NORMS OF THE COUPLING MATRICES 

WRITE (3,108) EST1,EST2,EST3,ESTA 

108 FORMAT ( 4E 16 . 6 ) 

EST=EST2*EST3 
ZET2=SMIN( 1 )*SMIN(2)/EST 
ZETA=SQRT(ZET2) 

WRITE(5, 105) 

WR I TE ( 5 , 308 ) ZETA 

308 FORMAT( ‘OTHE ESTIMATE OF THE COUPLING PARAMETAR IS ',E16.6) 

90 CALL EXIT 
END 


FEATURES SUPPORTED 
ONE WORD INTEGERS 
IOCS 

CORE REQUIREMENTS FOR 

COMMON 0 ... VARIABLE S 3 28 6 PROG RAM 1386 

END OF COMPILATION 
// XEQ 


66 



THE SYSTEM MATRIX IS 

0.0000 0.0463 0.438 7 0.0I3I 

0.8478 0.0000 0.0000 0.0000 

0.0000 0.0000 0.0000 1.0000 

0.0000 -1.0463 -3.2253 -0.0665 


THE ROOTS OF THE CHARACTERISTIC EQUATION ARE 

REAL PART IMAGINARY PART 

-0.000 5 0.290 9 

-0.0005 -0.2909 

-0.0327 1.7607 

-0.0327 -1.7607 

THE SIMILAR MATRIX IS 

-0.0005 0.2909 0.0000 0.0000 

r 0.2909 0.0005 0.0000 __0.0000_ 

0.0000 0.0000 -0.0327 1.7607 

0.0000 0.0000 -1.7607 -0.0327 


THE SIMILARITY TRANSFORMATION MATRIX IS 


0.3000 

0.8723 

-0.2888 

-0.0852 


0.3000 

-0.8757 

0.2935 

-0.0842 


0.1000 

0.0472 

-0.3919 

-0_._X267_ 


0.1000 
-0.0490 
0.4200 
M3. 7039 
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THE SYSTEM MATRIX IS 


0.0000 .. I . 0000 .... 0 . 0000 _ 0.0000 

-1.7865 -0.0534 0.0000 2.0000 

0.0000 0.0000 0.0000 1.0000 

0.0000 “2.0000 “904. 5314 -1.2036 


THE ROOTS OF THE CHARACTERISTIC EQUATION 


REAL PART IMAGINARY PART 

-0.0266 1.3334 

-0.0266 -1.3334 


-0.6019 30.1359 

-0.6019 -30.1359 

THE SIMILAR MATRIX IS 

-0.0266 1.3334 0.0000 0.0000 

-1.3334 ..-0,0266 0.0000 0.0000 

0.0000 0.0000 -0.6019 30.1359 

0.0000 0.0000 -30.1359 -0.6019 


THE SIMILARITY TRANSFORMATION MATRIX IS 

0.3000 0.3000 0.0030 

-0.4080 0.3920 -0.0939. ... 

0.0009 -0.0008 -0.0468 

0 . 001.1 0.0012 - L .3.18.6 


68 


ARE 


0.0030 

0.0894 

0.0446 

i.43_7J_ 



THE TRANSFORMED SELFCOUPLING MATRIX. _Q11 IS 

-0.3406 0.3652 -0.0 792 _0.0853 

-0.3632 0.3407 -0.0869 0.0772 

-0.1372 0.1210 -0.0332 0.0270 

-0.1084 0.1260 -0.0247 0.0299 


THE TRANSFORMED INTERCOUPLING MATRIX Q12 IS 

0.0495 
0.0488 
0.0178 
0 . 01 65 


THE TRANSFORMED. INTERCOUPLING MATRIX Q21 IS 

0.2307 
0.2291 
0.0738 
0.0865 


THE TRANSFORMED SELFCOUPLING MATRIX Q22 IS 


0.0354 • 0.0388 0.0000 0.0007 

-0.0369 -0.0374 -0.0002 -0.0004 

0.0717 -0.0674 0.0122 -0.0116 

0.0737 -0 . 0 72.1 0_. 0 1.28 -0.0122 


THE FOUR NORMS OF THE COUPLING MATRICES ARE 

0.766883E 00 0.631457E 00 0.145454E 01 0.144751E 00 


THE ESTIMATE OF THE COUPLING PARAMETAR IS 0.403822E-02 


0.6862 0 . 69 33 0 . 2 2.9 1. 

-0.6932 -0.6864 -0.2307 

0.3536 -0.3269 0.0848 

0.3431 -0.3698 0.0797 


-0.2826 
-0.3164 
-0.1219 
_— 0 .0870 


0.3063 
0.2690 
0.0925 
0 . 10_90 


-0.0517 
-0.0517. 
-0.0190 
-0 . 0_17JL_ 
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THE SYSTEM MATRIX IS 


0.0000 

0.0000 _ 

0 . 0000 _. 

1.0000 .... 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

1.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

1.0000 

0.0463 

- 1.2456 

- 0.4387 

- 1.2456 

1.0463 

- 0.0131 

0.0000 

- 0.8478 

0.0000 

- 0.1521 

0.0000 

0.0000 

1.0463 

- 1.2456 

- 3.2253 

- 1.2456 

1.0463 

- 0.0665 


THE ROOTS OF THE CHARACTERISTIC EQUATION ARE 
REAL PART I MAGI NAR Y~PART 

- 0. 1288 0.0000 

- 0. 8408 0.0000 

- 0.0783 1.720 4 

- 0.0783 - 1 . 720.4 

- 0.0929 0.8990 

- 0.0929 - 0.8990 

THE SIMILAR MATRIX IS 

- 0.1288 0.0000 0.0000 0.0000 0.0000 : 0.0000 

0.0000 - 0.8408 0.0000 0.0000 0...0000 0 . 0000 _ 

0.0000 0.0000 - 0.0783 1.7204 0.0000 0.0000 

0.0000 0.0000 r 1.7204 - 0.0783 0.0000 0.0000 

0.0000 0.0000 0.0000 0.0000 - 0.0929 0.8990 

0.0000 0.0000 0.0000 0.0000 - 0.8990 - 0 .. 0929 


THE SIMILARITY TRANSFORMATION MATRIX 


1.0000 1.0000 0.1000 

0 . 0 2 26 0.0822 - 0 . 0 11 2 .._ 

0.3635 0.4949 1.1220 

. - 0.1288 - 0.8408 - 0 . 1798 .. 

- 0.0029 - 0.0691 - 0.0220 

- 0.0468 - 0.4162 - 0.3200 .. 


IS 


0.1000 1.0000 1.0000 

_ 0 _. 01.33 0.9183 0.6 39 2 

0.1349 0.1836 0.0278 

0.1642 - 0.9919 . 0.8061 

- 0.0204 - 0.6600 0.7662 

1.9198 - 0.0420 .... 0.1625 
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THE SYSTEM MATRIX IS 


-0.1823 

0.0642 

0.0019.. 

0.0000 

0.0000 

0.0000 

0.0000 

1.0000 

0.0000 

0.0000 

0.1823 

-1.8508 

-0.0553 

0.0000 

2.0000 

0.0000 

. 0.0000 

0.0000 

0.0000 

1.0000 

-2.0000 

0.0000 

-2.0000 

-904.5314 

-1.2036 


THE ROOTS OF THE .CHARACTER I ST I.C EQUAT I ON. ARE 
RE ALP ART. IM AGIN A R Y PART 


-0.1761 

0.0000 




-0.0306 

1.3567 

-- - - . . 


-0.0306 

-1.3567 



-0.6019 

30.1359 



-0.6019 

-30.1359 


- - 

THE SIMILAR 

MATRIX IS 



-0. 1761 

0.0000 0.0000 

0.0000 

0.0000 

0.0000 

-0.0306 1.3567 

0.0000 

0.0000 

0.0000 

-1.3567 -0.0306 

0.0000 

0.0000 

0.0000 

0.0000 0.0000 

-0.6019 

30.1359 

0.0000 

0.0000 0.0000 

-30. 1359 

-0.6019 


SIMILARITY 

TRANSFORMATION 

MATRIX 

I S 


1 .0000 

0.0500 

0.0500 

0.0000 

0.0000 

0.0978 

-0.8765 

1.2306 

-0.0011 

0.0019 

-0.0172 

-1.6427 

1.2270 

-0.0588 

-0.0349 

-0.0021 

0.0035 ; 

0.0026 

-0.0293 

-0.0173 

0.0003 

-0.0036 

0.0047 

0.5413 

-0.8749 
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THE TRANSFORMED SELFCOUPLING MATRIX Qll IS 


-0.3595 

-0.1923 

-0.1968 

-0. 1009 

-0.1952 

-0.9725 

0.3918 

0.1893 

0.1559 

0.1159 

0.1537 

1 .1068 

-0.0990 

-0.0159 

-0.0163 

-0.0197 

-0.0159 

-0.1378 

-0.0580 

0.00-89 

-0.0156 

-0.0277 

-0.0196 

-0.2986 

-0 .0397 

-0.0996 

-0.0205 

-0.0007 

-0.0209 

-0.0217 

0.0126 

0.0533 

0.0196 

-0.0099 

0.0155 

-0.0739 


THE TRANSFORMED 

INTERCOUPLING MATRIX 

012 IS 


0.0189 

0.9356 

1.1327 

-9.0227 

-5.8359 

-0.0196 

-1.0610 

-1.2005 

9.9113 

6.9102 

0.0018 

0.1310 

0.1299 

-1.1356 

-0.7395 

0.0010 

0.2315 

0.1166 

-1.6091 

-1.0909 

0.0039 

0.0269 

0.1616 

-0.7517 

-Oi 9860 

-0.0030 

0.0621 

-0.1180 

0.1376 

0.0888 


THE TRANSFORMED 

INTERCOUPLING MATRIX 

021 IS 



-0.0068 

-0.0732 

-0.0281 

0.0269 

0.2685 

0.3809 

0.0793 

0.5719 

0.1935 

-0.1338 

0.0695 

-1.0890 

0.0588 

0.9139 

0.1050 

. -0.0968 

0.0251 

-0.7966 

0.9999 

0.1576 

0.1815 

0.1691 

0.1979 

1.5636 

-0.8110 

-0.0792 

-0.2589 _ 

-0.3311 

-0.2659 

-3.0315 


THE TRANSFORMED 

SELFCOUPLING 

MATRIX Q 22 

IS 


0.0000 

-0.0000 

0.0001 . 

-0.0000 

-0.0000 

0.0019 

-0.0373 

0.0507 

-0.0053 

-0.0033 

0.0009 

-0.0280 

0.0350 

0.0076 

0.0050 

-0.0037 

-0.9056 

•0.2929 

3.1593 

2.0903 

0.0067 

... 0 . 6 J..9 0 . 

_0„.3 1 00 

-5.3337 

-3.9999 


THE 

FOUR NORMS 

OF THE COUPLING MATRICES ARE 


0 

. 165203E 01 

0.192159E 02 0.379971E 

01 0 • 799500E 01 





THE 

ESTIMATE OF 

THE COUPLING PARAMETAR IS 

0.671987E-02 


NASA-Langley, 1974 
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